/*
 * topological-hex, a program to compute combinatorial hexahedral meshes.
 *
 * Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
 *
 * List of the contributors to the development, description and complete
 * License: see LICENSE file.
 *
 * This program (topological-hex) is free software:
 * you can redistribute it and/or modify it under the terms
 * of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program (see COPYING file).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "dfs.h"
#include "shellable_dfs.h"
#include "symmetry.h"
#include "lower_bound.h"
#include "homology.h"

/**
 * Explores a sub-tree, where @p id is the number of hexahedra that have
 * been completed so far.
 */
static void start_hex(solver *solver);

/**
 * Restricts te domain of variable @p var to the vertices contained in @p
 * max_domain.
 */
static void restrict_domain(solver *solver, fixed_point_delta *delta,
                            quad_vertex_index var, const vertex_set *max_domain);

/**
 * Hook called whenever both vertices of an edge have been chosen.
 *
 * This updates the set of allowed neighbors in order to avoid creating edges
 * which correspond to diagonals of other hexahedra. If the created edge is (@p
 * a, @p b), then any edge (@p a, v) where v is a neighbor of @p b, is
 * forbidden.
 *
 * @param i Index of the first vertex within the hexahedron
 * @param j Index of the second vertex within the hexahedron
 * @param a Index of the first vertex
 * @param b Index of the second vertex
 */
static void finish_edge(solver *solver, fixed_point_delta *delta,
                        hex_vertex_index i, hex_vertex_index j,
                        vertex_index a, vertex_index b);

/**
 * Guarantees that the edge (@p a, @p b) can no longer be created.
 *
 * This is called for pairs of non-adjacent vertices within a hexahedron.
 */
static void forbid_edge(solver *solver, fixed_point_delta *delta,
                        hex_vertex_index i, hex_vertex_index j,
                        vertex_index a, vertex_index b);

/**
 * Checks whether or not a pair of vertices can be used as the diagonal of a
 * hexahedron or of a facet.
 *
 * This is true if (@p a, @p b) is neither a diagonal nor an edge of any
 * existing hexahedron.
 */
static bool check_cube_diagonal(solver *solver, fixed_point_delta *delta,
                                hex_index hex_id,
                                vertex_index a, vertex_index b);

/**
 * Hook called whenever a diagonal of a hexahedron is finished.
 *
 * This checks that creating this diagonal is allowed.
 *
 * @param i Index of the vertex within the top face of the hexahedron.
 * @param a First vertex of the diagonal
 * @param b Second vertex of the diagonal
 *
 * @return True if adding this vertex does not create any incompatibilities
 */
static bool finish_cube_diagonal(solver *solver, fixed_point_delta *delta,
                                 hex_index hex_id, quad_vertex_index i,
                                 vertex_index a, vertex_index b);

/**
 * Hook called whenever the diagonal of one of the facets on the side of the
 * hexahedron has been completed.
 *
 * If this diagonal matches an existing edge, the search has failed. If it
 * matches the diagonal of another face, this can be used to determine the value
 * of variable @p other_var.
 *
 * @param a First vertex of the diagonal
 * @param b First vertex of the diagonal
 * @param other_var Index of the variable on the same facet whose value is not
 *   necessarily yet bound. There will always be the two vertices on the
 *   diagonal, one vertex on the bottom of the hexahedron, and a 4th variable
 *   whose domain is going to be restricted by this function call.
 *
 * @return True if creating this facet was allowed.
 */
static bool finish_diagonal(solver *solver, fixed_point_delta *delta,
                            hex_index hex_id,
                            vertex_index a, vertex_index b,
                            quad_vertex_index other_var);

/**
 * Hook called whenever the diagonal of the top facet has been completed.
 *
 * This behaves the same as @ref finish_diagonal, except there may be two
 * vertices whose values are not currently known.
 */
static bool finish_top_diagonal(
  solver *solver, fixed_point_delta *delta,
  hex_index hex_id,
  vertex_index a, vertex_index b,
  quad_vertex_index var_a, quad_vertex_index var_b);

/**
 * Returns a set containing the two vertices on the diagonal of a quad described
 * by @p descriptor.
 *
 * A quad descriptor only encodes one of the diagonals, the other diagonal is
 * used as an index when looking it up.
 */
static void quad_descriptor_diagonal_set(quad_descriptor descriptor,
                                         vertex_set *set);

/**
 * Returns a set containing the four vertices opposite to a quad if it has
 * already appeared in a hexahedron. Otherwise, returns an empty set.
 */
static void quad_opposite_set(solver *solver, quad_descriptor desc,
                              vertex_set *set);

/**
 * Checks whether or not two quadrangles (which this function assumes contain
 * the same 4 vertices) have the same orientation.
 */
static bool quad_has_same_orientation(vertex_index *q_a, vertex_index *q_b);

/**
 * Hook called whenever a quad has been completed.
 *
 * The vertices of the quad are (@p a, @p c, @p b, @p d). In other words, the
 * first diagonal is (@p a, @p b) and the second diagonal is (@p c, @p d).
 *
 * @return True if creating this quad was allowed.
 */
static bool finish_quad(solver *solver, fixed_point_delta *delta,
                        hex_index hex_id,
                        vertex_index a, vertex_index b,
                        vertex_index c, vertex_index d);

void solver_run(solver *solver, solver_callback cb, void *user_data) {
  solver->cb        = cb;
  solver->user_data = user_data;

  if (solver->flags & SOLVER_SHELLABLE_ONLY)
    shellable_dfs(solver);
  else
    start_hex(solver);
}

/**
 * Removes an element from the quad queue.
 *
 * This is used to cancel the appending of quads when backtracking from the node
 * where they first occurred.
 */
static void quad_queue_remove(solver *solver, quad_index q);

/**
 * Initializes a new hexahedron whose bottom face is @p quad.
 *
 * This function computes the initial domain for the 4 vertices whose values are
 * not determined by @p quad.
 */
static void init_domains(solver *solver, hex_index hex_id,
                         quad_index quad);

static void start_hex(solver *solver) {
  solver->hex_id++;

  uint32_t delta_queue;
  quad_index quad = pop_quad(solver, &delta_queue);

  if (quad == NO_QUAD) {
    if (solver->quad_queue_offset == solver->quad_queue_size)
      solver->cb(solver, solver->hex_id, solver->user_data);
  }
  else if (solver->hex_id != solver->num_hexes) {
    init_domains(solver, solver->hex_id, quad);

    fixed_point_delta trail;
    if (fixed_point(solver, &trail)) {
#ifdef VERBOSE
      char header[256];
      snprintf(header, sizeof(header), "Starting hex #%u", solver->hex_id);
      solver_print(solver, solver->hex_id, header);
#endif

      choose_vertex(solver);
    }
    undo_fixed_point(solver, &trail);
  }

  undo_pop_quad(solver, quad, delta_queue);
  solver->hex_id--;
}

void choose_vertex(solver *solver) {
  vertex_set restricted_domain;
  quad_vertex_index vertex_id = choose_variable(solver, &restricted_domain);

  if (vertex_id == 4) {
    ALIGNED vertex_set domains[4];
    uint32_t bound_mask = solver->bound_mask;
    for (quad_vertex_index i = 0; i < 4; i++) domains[i] = solver->domains[i];
    start_hex(solver);
    for (quad_vertex_index i = 0; i < 4; i++) solver->domains[i] = domains[i];
    solver->bound_mask = bound_mask;
  }
  else {
    vertex_set old_domain = solver->domains[vertex_id];

    solver->num_branches++;

#ifdef VERBOSE
    printf("domain (%u, %u): {", solver->hex_id, vertex_id);
    VERTEX_SET_FOREACH(restricted_domain, v) {
      printf("%u, ", (vertex_index)v);
    }
    printf("}\n");
#endif

    no_goods_delta no_goods_trail;
    solver_start_branch(solver, &no_goods_trail);

    VERTEX_SET_FOREACH(restricted_domain, v) {
      vertex_set_init_empty(solver->domains + vertex_id);
      vertex_set_insert(solver->domains + vertex_id, v);

      bool increased = false;
      if (v >= solver->max_vertex_id) {
        solver->max_vertex_id++;
        increased = true;
      }

#ifdef VERBOSE
      char pre_header[256];
      snprintf(pre_header, sizeof(pre_header),
               "Before choosing Hex #%u, vertex #%u = %u",
               solver->hex_id, vertex_id, (vertex_index)v);
      solver_print(solver, solver->hex_id, pre_header);
#endif

      fixed_point_delta trail;
      if (fixed_point(solver, &trail)) {
#ifdef VERBOSE
        char header[256];
        snprintf(header, sizeof(header),
                 "After choosing Hex #%u, vertex #%u = %u",
                 solver->hex_id, vertex_id, (vertex_index)v);
        solver_print(solver, solver->hex_id, header);
#endif
        choose_vertex(solver);
      }
      else {
#ifdef VERBOSE
        char header[256];
        snprintf(header, sizeof(header),
                 "After failing Hex #%u, vertex #%u = %u",
                 solver->hex_id, vertex_id, (vertex_index)v);
        solver_print(solver, solver->hex_id, header);
#endif
      }

      undo_fixed_point(solver, &trail);
      if (increased) solver->max_vertex_id--;

      solver_is_no_good(solver, vertex_id, v, &no_goods_trail);

#ifdef VERBOSE
      char header[256];
      snprintf(header, sizeof(header),
               "After backtrack from Hex #%u, vertex #%u = %u",
               solver->hex_id, vertex_id, (vertex_index)v);
      solver_print(solver, solver->hex_id, header);
#endif
    }

    solver_undo_no_goods(solver, &no_goods_trail);
    solver->domains[vertex_id] = old_domain;
  }
}

bool propagate(solver *solver) {
  bool success = true;

  fixed_point_delta delta;

  while (true) {
    vertex_set domain;
    quad_vertex_index vertex_id = choose_variable(solver, &domain);
    if (vertex_id == 4) {
      solver->hex_id++;

      uint32_t num_quads;
      quad_index quad = pop_quad(solver, &num_quads);
      if (quad == NO_QUAD) {
        if (solver->quad_queue_offset == solver->quad_queue_size)
          solver->cb(solver, solver->hex_id, solver->user_data);
        success = false;
        break;
      }
      else if (solver->hex_id == solver->num_hexes) {
        success = false;
        break;
      }

      init_domains(solver, solver->hex_id, quad);

      if (!fixed_point(solver, &delta)) {
        success = false;
        break;
      }

      vertex_id = 0;
    }
    else {
      if (solver->bound_mask & (1 << vertex_id))
        vertex_id++;
      else {
        vertex_set old_domain = solver->domains[vertex_id];
        solver->domains[vertex_id] = domain;

        if (!fixed_point(solver, &delta)) {
          success = false;
          solver->domains[vertex_id] = old_domain;
          break;
        }

        if (vertex_set_unit(&domain)) {
          if (vertex_set_min(&domain) >= solver->max_vertex_id)
            solver->max_vertex_id++;
          vertex_id++;
        }
        else
          break;
      }
    }
  }

  return success;
}

bool fixed_point(solver *solver, fixed_point_delta *delta) {
  bool failed = false;

  uint32_t bound_mask = solver->bound_mask;

  hex_index hex_id = solver->hex_id;

  for (quad_vertex_index i = 0; i < 4; i++) {
    if ((solver->bound_mask & (1 << i)) == 0 &&
        vertex_set_unit(&solver->domains[i])) {
      bound_mask |= (1 << i);
      solver->hexes[8*hex_id + 4 + i] = vertex_set_min(solver->domains + i);
    }
    else if (vertex_set_empty(solver->domains + i))
      failed = true;
  }

  delta->new_bound_mask = bound_mask;
  delta->diagonal_mask = 0;
  for (quad_vertex_index i = 0; i < 4; i++)
    vertex_set_init_empty(delta->domain_changes + i);

  delta->max_quad_change = 0;
  delta->num_occurrences_changes = 0;

  for (hex_vertex_index i = 0; i < 8; i++)
    vertex_set_init_empty(delta->removed_from + i);
  for (hex_vertex_index i = 0; i < 8; i++)
    vertex_set_init_empty(delta->allowed_set_changes + i);
  for (hex_vertex_index i = 0; i < 8; i++)
    vertex_set_init_empty(delta->new_neighbors + i);

  vertex_set_init_empty(&delta->new_even_vertices);
  vertex_set_init_empty(&delta->new_odd_vertices);

  uint32_t to_visit = bound_mask & ~solver->bound_mask;
  while (!failed && to_visit != 0) {
    quad_vertex_index i = find_lsb32(to_visit);
    to_visit &= ~(1 << i);

    quad_vertex_index succ_i = (i + 1) % 4;
    quad_vertex_index pred_i = (i + 3) % 4;
    quad_vertex_index opposite_i = (i + 2) % 4;

    vertex_index v = solver->hexes[8*hex_id + 4+i];
    vertex_index below_v = solver->hexes[8*hex_id + i];

    vertex_index succ_v = solver->hexes[8*hex_id + 4+succ_i];
    vertex_index below_succ_v = solver->hexes[8*hex_id + succ_i];

    vertex_index pred_v = solver->hexes[8*hex_id + 4+pred_i];
    vertex_index below_pred_v = solver->hexes[8*hex_id + pred_i];

    vertex_index opposite_v = solver->hexes[8*hex_id + 4+opposite_i];
    vertex_index below_opposite_v = solver->hexes[8*hex_id + opposite_i];

    if (!finish_cube_diagonal(solver, delta, hex_id, i, v, below_opposite_v)) {
      failed = true;
      break;
    }

    finish_edge(solver, delta, i, 4 + i, below_v, v);

    if (bound_mask & (1 << succ_i))
      finish_edge(solver, delta, 4 + succ_i, 4 + i, succ_v, v);

    if (!finish_diagonal(solver, delta, hex_id, v, below_succ_v, succ_i)) {
      failed = true;
      break;
    }

    if (bound_mask & (1 << pred_i))
      finish_edge(solver, delta, 4 + pred_i, 4 + i, pred_v, v);

    if (!finish_diagonal(solver, delta, hex_id, v, below_pred_v, pred_i)) {
      failed = true;
      break;
    }

    forbid_edge(solver, delta, 4+i, opposite_i, v, below_opposite_v);
    forbid_edge(solver, delta, 4+i, pred_i, v, below_pred_v);
    forbid_edge(solver, delta, 4+i, succ_i, v, below_succ_v);

    if (solver->quad_descriptors[
          diagonal_index(solver->num_vertices, v, below_opposite_v)] !=
        EMPTY_QUAD_DESCRIPTOR) {
      failed = true;
      break;
    }

    if (bound_mask & (1 << opposite_i)) {
      forbid_edge(solver, delta, 4 + i, 4 + opposite_i, v, opposite_v);
      if (!finish_top_diagonal(solver, delta, hex_id, v, opposite_v,
                               succ_i, pred_i)) {
        failed = true;
        break;
      }
    }

    for (quad_vertex_index j = 0; j < 4; j++) {
      if (j != i) {
        if (vertex_set_remove(solver->domains + j, v))
          vertex_set_insert(delta->domain_changes + j, v);
      }

      quad_vertex_index succ_j = (j+1) % 4;
      quad_vertex_index pred_j = (j+3) % 4;
      quad_vertex_index opposite_j = (j+2) % 4;

      vertex_set to_keep = solver->domains[j];
      if ((bound_mask & (1 << succ_j)) != 0) {
        vertex_set_intersect(&to_keep,
                             solver->allowed_neighbors +
                             solver->hexes[8*hex_id + 4+succ_j]);
      }

      if ((bound_mask & (1 << pred_j)) != 0) {
        vertex_set_intersect(&to_keep,
                             solver->allowed_neighbors +
                             solver->hexes[8*hex_id + 4+pred_j]);
      }

      if ((bound_mask & (1 << opposite_j)) != 0) {
        vertex_set_subtract(&to_keep,
                            solver->known_diagonals +
                            solver->hexes[8*hex_id + 4+opposite_j]);
      }

      vertex_set_intersect(&to_keep,
                           solver->allowed_neighbors +
                           solver->hexes[8*hex_id + j]);

      vertex_set removed = solver->domains[j];
      vertex_set_subtract(&removed, &to_keep);

      vertex_set_union(delta->domain_changes + j, &removed);
      solver->domains[j] = to_keep;

      if ((bound_mask & (1 << j)) == 0 &&
          vertex_set_unit(solver->domains + j)) {
        bound_mask |= (1 << j);
        to_visit |= (1 << j);
        solver->hexes[8*hex_id + 4+j] = vertex_set_min(solver->domains + j);
      }
      else if (vertex_set_empty(solver->domains + j)) {
        failed = true;
        break;
      }
    }

   if (failed) break;

    if ((bound_mask & (1 << succ_i)) != 0 &&
        (bound_mask & (1 << (4+i))) == 0) {
      bound_mask |= (1 << (4+i));

      vertex_index a = solver->hexes[8*hex_id + 4+i];
      vertex_index b = solver->hexes[8*hex_id + succ_i];
      vertex_index c = solver->hexes[8*hex_id + i];
      vertex_index d = solver->hexes[8*hex_id + 4+succ_i];

      if (!finish_quad(solver, delta, hex_id, a, b, c, d)) {
        failed = true;
        break;
      }
    }

    if ((bound_mask & (1 << pred_i)) != 0 &&
        (bound_mask & (1 << (4+pred_i))) == 0) {
      bound_mask |= (1 << (4+((i+3)%4)));

      vertex_index a = solver->hexes[8*hex_id + 4+i];
      vertex_index b = solver->hexes[8*hex_id + pred_i];
      vertex_index c = solver->hexes[8*hex_id + 4+pred_i];
      vertex_index d = solver->hexes[8*hex_id + i];

      if (!finish_quad(solver, delta, hex_id, a, b, c, d)) {
        failed = true;
        break;
      }
    }

    if ((bound_mask & ((1 << 4)-1)) == (1 << 4)-1 &&
        (bound_mask & (1 << 8)) == 0) {
      bound_mask |= (1 << 8);

      vertex_index a = solver->hexes[8*hex_id + 4];
      vertex_index b = solver->hexes[8*hex_id + 6];
      vertex_index c = solver->hexes[8*hex_id + 5];
      vertex_index d = solver->hexes[8*hex_id + 7];

      if (!finish_quad(solver, delta, hex_id, a, b, c, d)) {
        failed = true;
        break;
      }
    }
  }

  delta->new_bound_mask = bound_mask & ~solver->bound_mask;
  solver->bound_mask = bound_mask;

  for (quad_vertex_index i = 0; i < 4; i++) {
    if (delta->new_bound_mask & (1 << i)) {
      vertex_index v = solver->hexes[8 * hex_id + 4 + i];
      if (vertex_set_contains(&solver->even_vertices,
                              solver->hexes[8 * hex_id + i])) {
        if (!vertex_set_contains(&solver->odd_vertices, v)) {
          vertex_set_insert(&delta->new_odd_vertices, v);
          vertex_set_insert(&solver->odd_vertices, v);
        }
      }
      else {
        if (!vertex_set_contains(&solver->even_vertices, v)) {
          vertex_set_insert(&solver->even_vertices, v);
          vertex_set_insert(&delta->new_even_vertices, v);
        }
      }
    }
  }

  if (!failed) {
    uint32_t next_quad;
    for (next_quad = solver->quad_queue_offset;
         next_quad < solver->quad_queue_size &&
           solver->num_occurrences[solver->quad_queue[next_quad]] == 0;
         next_quad++);

    if (next_quad == solver->quad_queue_size ||
        solver->hex_id % 8 == 0)
      failed = delta->new_bound_mask != 0 && is_dominated(solver);
  }

  if (!failed && (solver->bound_mask & ((1 << 4)-1)) == (1 << 4) - 1)
    failed = !check_homology(solver);

  return !failed;
}

void undo_fixed_point(solver *solver, const fixed_point_delta *delta) {
  hex_index hex_id = solver->hex_id;

  solver->bound_mask &= ~delta->new_bound_mask;

  for (uint32_t i = 0; i < delta->num_occurrences_changes; i++) {
    quad_index q = delta->occurred[i];
    solver->num_occurrences[q]++;
    if (solver->num_occurrences[q] == 2)
      quad_queue_remove(solver, q);
  }

  for (quad_index i = solver->max_quad_id - delta->max_quad_change;
       i < solver->max_quad_id; i++) {
    vertex_index a = solver->quads[4*i + 0];
    vertex_index b = solver->quads[4*i + 1];
    vertex_index c = solver->quads[4*i + 2];
    vertex_index d = solver->quads[4*i + 3];

    solver->quad_descriptors[diagonal_index(solver->num_vertices, a, c)] =
      EMPTY_QUAD_DESCRIPTOR;
    solver->quad_descriptors[diagonal_index(solver->num_vertices, b, d)] =
      EMPTY_QUAD_DESCRIPTOR;

    vertex_set_remove(solver->known_quad_diagonals + a, c);
    vertex_set_remove(solver->known_quad_diagonals + c, a);

    vertex_set_remove(solver->known_quad_diagonals + b, d);
    vertex_set_remove(solver->known_quad_diagonals + d, b);
  }
  solver->max_quad_id -= delta->max_quad_change;

  for (hex_vertex_index i = 0; i < 8; i++) {
    if (i >= 4) {
      vertex_set_union(solver->domains + (i - 4),
                       delta->domain_changes + (i - 4));

    }

    VERTEX_SET_FOREACH(delta->removed_from[i], x) {
      vertex_set_insert(solver->allowed_neighbors + x,
                        solver->hexes[8*hex_id + i]);
    }

    if (!vertex_set_empty(delta->allowed_set_changes + i)) {
      vertex_set_union(solver->allowed_neighbors +
                       solver->hexes[8*hex_id + i],
                       delta->allowed_set_changes + i);
    }

    if (!vertex_set_empty(delta->new_neighbors + i)) {
      vertex_set_subtract(solver->known_neighbors +
                          solver->hexes[8*hex_id + i],
                          delta->new_neighbors + i);
    }
  }

  for (quad_vertex_index i = 0; i < 4; i++) {
    if (delta->diagonal_mask & (1 << i)) {
      vertex_index v = solver->hexes[8*hex_id + 4 + i];
      vertex_index u = solver->hexes[8*hex_id + (i+2)%4];
      vertex_set_remove(solver->known_diagonals + v, u);
      vertex_set_remove(solver->known_diagonals + u, v);
    }
  }

  vertex_set_subtract(&solver->even_vertices, &delta->new_even_vertices);
  vertex_set_subtract(&solver->odd_vertices,  &delta->new_odd_vertices);
}

static
void recompute_finished_vertices(solver *solver) {
  vertex_set_init_up_to(&solver->finished_vertices,
                        solver->max_vertex_id);

  for (uint32_t i = solver->quad_queue_offset;
       i < solver->quad_queue_size; i++) {
    quad_index q = solver->quad_queue[i];
    if (solver->num_occurrences[q] == 0) continue;

    for (quad_vertex_index j = 0; j < 4; j++)
      vertex_set_remove(&solver->finished_vertices, solver->quads[4 * q + j]);
  }
}

static
bool compute_locally_finished_vertices(solver *solver, quad_index q) {
  vertex_set_init_up_to(&solver->finished_vertices,
                        solver->max_vertex_id);

  vertex_set seen;
  vertex_set_init_empty(&seen);
  for (quad_vertex_index i = 0; i < 4; i++)
    vertex_set_insert(&seen, solver->quads[4 * q + i]);

  uint32_t offset = solver->quad_queue_offset + 1;
  while (true) {
    uint32_t i;
    for (i = offset; i < solver->quad_queue_size; i++) {
      quad_index q = solver->quad_queue[i];
      if (solver->num_occurrences[q] == 0) continue;

      quad_vertex_index j;
      for (j = 0; j < 4 &&
             !vertex_set_contains(&seen, solver->quads[4 * q + j]);
           j++);
      if (j != 4) {
        break;
      }
    }

    if (i == solver->quad_queue_size) break;

    quad_index q = solver->quad_queue[i];
    quad_index to_move = solver->quad_queue[offset];

    solver->quad_queue[offset++] = q;
    solver->quad_queue[i] = to_move;

    for (quad_vertex_index j = 0; j < 4; j++)
      vertex_set_insert(&seen, solver->quads[4 * q + j]);
  }

  vertex_set_subtract(&solver->finished_vertices, &seen);

  uint32_t num_quads = offset - solver->quad_queue_offset;
  uint32_t num_vertices = solver->num_vertices - solver->max_vertex_id;
  uint32_t num_hexes = solver->num_hexes - solver->hex_id;

  vertex_index min_v = 0;
  hex_index min_h = 0;

  component_lower_bound(
    solver, solver->quad_queue + solver->quad_queue_offset,
    num_quads, &min_v, &min_h);

  if (num_vertices < min_v || num_hexes < min_h)
    return false;

  return true;
}

quad_index pop_quad(solver *solver, uint32_t *delta_size) {
  recompute_finished_vertices(solver);

  uint32_t old_offset = solver->quad_queue_offset;
  while (solver->quad_queue_offset < solver->quad_queue_size &&
         solver->num_occurrences[
           solver->quad_queue[solver->quad_queue_offset]] == 0) {
    solver->quad_queue_offset++;
  }

  *delta_size = solver->quad_queue_offset - old_offset;
  if (solver->quad_queue_offset == solver->quad_queue_size)
    return NO_QUAD;

  size_t max_size = SIZE_MAX;
  uint32_t best_i = UINT32_MAX;

  vertex_set even = solver->even_vertices;
  vertex_set odd = solver->odd_vertices;

  vertex_set_init_empty(&solver->even_vertices);
  vertex_set_init_empty(&solver->odd_vertices);

  for (uint32_t delta = 0;
       delta < solver->quad_queue_size - solver->quad_queue_offset; delta++) {
    /* Consider the most recently added quads first seems to be faster */
    uint32_t i = solver->quad_queue_size - delta - 1;
    /* uint32_t i = solver->quad_queue_offset + delta; */

    quad_index q = solver->quad_queue[i];
    if (solver->num_occurrences[q] == 0) continue;

    init_domains(solver, solver->hex_id, q);

    fixed_point_delta delta;

    bool do_fix_point = false;
    for (quad_vertex_index j = 0; j < 4; j++)
      do_fix_point |= vertex_set_unit(solver->domains + j);

    if (do_fix_point) {
      if (!fixed_point(solver, &delta)) {
        best_i = i;
        max_size = 0;
        undo_fixed_point(solver, &delta);
        break;
      }
    }

    size_t size = 1;
    for (quad_vertex_index j = 0; j < 4; j++) {
      vertex_set small_domain;
      vertex_set_init_up_to(
        &small_domain, solver->max_vertex_id == solver->num_vertices ?
        solver->max_vertex_id : solver->max_vertex_id + 1);
      vertex_set_intersect(&small_domain, solver->domains + j);

      size *= vertex_set_size(&small_domain);
    }

    if (do_fix_point)
      undo_fixed_point(solver, &delta);

    if (size < max_size) {
      max_size = size;
      best_i = i;
      if (max_size <= 1) break;
    }
  }

  solver->even_vertices = even;
  solver->odd_vertices  = odd;

  quad_index old_first = solver->quad_queue[solver->quad_queue_offset];
  quad_index best_quad = solver->quad_queue[best_i];

  solver->quad_queue[solver->quad_queue_offset] = best_quad;
  solver->quad_queue[best_i] = old_first;

  if (!compute_locally_finished_vertices(solver, best_quad))
    return NO_QUAD;

  solver->quad_queue_offset++;
  (*delta_size)++;
  solver->num_occurrences[best_quad]--;

  return best_quad;
}

void undo_pop_quad(solver *solver, quad_index quad,
                          uint32_t delta_size) {
  solver->quad_queue_offset -= delta_size;
  if (quad != NO_QUAD)
    solver->num_occurrences[quad]++;
}

static void quad_queue_remove(solver *solver, quad_index q) {
  uint32_t i = 0;
  for (i = solver->quad_queue_offset;
       i < solver->quad_queue_size && solver->quad_queue[i] != q; i++);

  solver->quad_queue[i] = solver->quad_queue[solver->quad_queue_size - 1];
  solver->quad_queue_size--;
}

static void init_domains(solver *solver, hex_index hex_id, quad_index quad) {
  memcpy(solver->hexes + 8*hex_id, solver->quads + 4*quad,
         4*sizeof(*solver->hexes));

  vertex_index tmp = solver->hexes[8*hex_id + 1];
  solver->hexes[8*hex_id + 1] = solver->hexes[8*hex_id + 3];
  solver->hexes[8*hex_id + 3] = tmp;

  vertex_set opposite;
  quad_opposite_set(
    solver,
    solver->quad_descriptors[
      diagonal_index(
        solver->num_vertices,
        solver->quads[4*quad + 0], solver->quads[4*quad + 2])],
    &opposite);

  for (quad_vertex_index i = 0; i < 4; i++) {
    solver->domains[i] =
      solver->allowed_neighbors[solver->hexes[8*hex_id + i]];
    vertex_set_subtract(solver->domains + i, &solver->finished_vertices);

    if (solver->flags & SOLVER_BIPARTITE_GRAPH) {
      if (vertex_set_contains(&solver->even_vertices,
                              solver->hexes[8 * hex_id + i]))
        vertex_set_subtract(solver->domains + i, &solver->even_vertices);
      else
        vertex_set_subtract(solver->domains + i, &solver->odd_vertices);
    }

    for (quad_vertex_index j = 0; j < 4; j++) {
      vertex_set_remove(solver->domains + i, solver->hexes[8*hex_id + j]);

      if (i != j) {
        vertex_set_subtract(solver->domains + i, solver->known_neighbors +
                            solver->hexes[8*hex_id + j]);
        vertex_set_subtract(solver->domains + i, solver->known_diagonals +
                            solver->hexes[8*hex_id + j]);
      }
    }

    vertex_set_subtract(solver->domains + i, &opposite);
    vertex_set_subtract(
      solver->domains + i,
      solver->known_quad_diagonals + solver->hexes[8*hex_id + (i+2)%4]);
  }

  solver->bound_mask = 0;
}

quad_vertex_index choose_variable(const solver *solver, vertex_set *domain) {
  quad_vertex_index best_id = 4;
  vertex_index best_size = (vertex_index)-1;

  for (quad_vertex_index i = 0; i < 4; i++) {
    if (solver->bound_mask & (1 << i)) continue;

    vertex_set restricted_domain = solver->domains[i];
    vertex_set tmp;
    vertex_set_init_up_to(
      &tmp, solver->max_vertex_id == solver->num_vertices ?
      solver->max_vertex_id : solver->max_vertex_id + 1);
    vertex_set_intersect(&restricted_domain, &tmp);

    vertex_index size = vertex_set_size(&restricted_domain);
    if (size < best_size) {
      *domain = restricted_domain;
      best_id = i;
      best_size = size;
    }
  }

  return best_id;
}

static void restrict_domain(solver *solver, fixed_point_delta *delta,
                            quad_vertex_index var,
                            const vertex_set *max_domain) {
  vertex_set removed = solver->domains[var];
  vertex_set_subtract(&removed, max_domain);
  vertex_set_union(delta->domain_changes + var, &removed);
  vertex_set_intersect(solver->domains + var, max_domain);
}

static void finish_ordered_edge(solver *solver, fixed_point_delta *delta,
                                hex_vertex_index i, hex_vertex_index j,
                                vertex_index a, vertex_index b);

static void finish_edge(solver *solver, fixed_point_delta *delta,
                        hex_vertex_index i, hex_vertex_index j,
                        vertex_index a, vertex_index b) {
  finish_ordered_edge(solver, delta, i, j, a, b);
  finish_ordered_edge(solver, delta, j, i, b, a);
}

static void finish_ordered_edge(solver *solver, fixed_point_delta *delta,
                                hex_vertex_index i, hex_vertex_index j,
                                vertex_index a, vertex_index b) {
  if (vertex_set_insert(solver->known_neighbors + a, b)) {
    vertex_set_insert(delta->new_neighbors + i, b);
    VERTEX_SET_FOREACH(solver->known_neighbors[a], v) {
      if (vertex_set_remove(solver->allowed_neighbors + v, b))
        vertex_set_insert(delta->removed_from + j, v);
    }

    vertex_set removed = solver->allowed_neighbors[a];
    vertex_set_intersect(&removed, solver->known_neighbors + b);
    vertex_set_union(delta->allowed_set_changes + i, &removed);
    vertex_set_subtract(solver->allowed_neighbors + a, &removed);
  }
}

static void forbid_ordered_edge(solver *solver, fixed_point_delta *delta,
                                hex_vertex_index i, hex_vertex_index j,
                                vertex_index a, vertex_index b);

static void forbid_edge(solver *solver, fixed_point_delta *delta,
                        hex_vertex_index i, hex_vertex_index j,
                        vertex_index a, vertex_index b) {
  forbid_ordered_edge(solver, delta, i, j, a, b);
  forbid_ordered_edge(solver, delta, j, i, b, a);
}

static void forbid_ordered_edge(solver *solver, fixed_point_delta *delta,
                                hex_vertex_index i, hex_vertex_index j,
                                vertex_index a, vertex_index b) {
  if (vertex_set_remove(solver->allowed_neighbors + a, b))
    vertex_set_insert(delta->removed_from + i, b);
}

static bool check_cube_diagonal(solver *solver, fixed_point_delta *delta,
                                hex_index hex_id,
                                vertex_index a, vertex_index b) {
  return !vertex_set_contains(solver->known_diagonals + a, b) &&
    !vertex_set_contains(solver->known_neighbors + a, b);
}

static bool finish_cube_diagonal(solver *solver, fixed_point_delta *delta,
                                 hex_index hex_id, quad_vertex_index i,
                                 vertex_index a, vertex_index b) {
  if (!check_cube_diagonal(solver, delta, hex_id, a, b) ||
      solver->quad_descriptors[diagonal_index(solver->num_vertices, a, b)] !=
      EMPTY_QUAD_DESCRIPTOR) {
    return false;
  }
  else {
    delta->diagonal_mask |= (1 << i);
    vertex_set_insert(solver->known_diagonals + a, b);
    vertex_set_insert(solver->known_diagonals + b, a);
  }

  return true;
}

static bool finish_diagonal(solver *solver, fixed_point_delta *delta,
                            hex_index hex_id,
                            vertex_index a, vertex_index b,
                            quad_vertex_index other_var) {
  if (!check_cube_diagonal(solver, delta, hex_id, a, b))
    return false;

  quad_descriptor quad = solver->quad_descriptors[
    diagonal_index(solver->num_vertices, a, b)];

  if (quad != EMPTY_QUAD_DESCRIPTOR) {
    vertex_set dom;
    quad_descriptor_diagonal_set(quad, &dom);
    restrict_domain(solver, delta, other_var, &dom);
  }

  return true;
}

static bool finish_top_diagonal(
  solver *solver, fixed_point_delta *delta,
  hex_index hex_id,
  vertex_index a, vertex_index b,
  quad_vertex_index var_a, quad_vertex_index var_b) {
  if (!check_cube_diagonal(solver, delta, hex_id, a, b))
    return false;

  quad_descriptor quad = solver->quad_descriptors[
    diagonal_index(solver->num_vertices, a, b)];

  if (quad != EMPTY_QUAD_DESCRIPTOR) {
    vertex_set dom;
    quad_descriptor_diagonal_set(quad, &dom);
    restrict_domain(solver, delta, var_a, &dom);
    restrict_domain(solver, delta, var_b, &dom);
  }

  return true;
}

static void quad_descriptor_diagonal_set(quad_descriptor descriptor,
                                         vertex_set *set) {
  vertex_set_init_empty(set);

  if (quad_descriptor_first_vertex(descriptor) != VERTEX_SET_MAX)
    vertex_set_insert(set, quad_descriptor_first_vertex(descriptor));

  if (quad_descriptor_second_vertex(descriptor) != VERTEX_SET_MAX)
    vertex_set_insert(set, quad_descriptor_second_vertex(descriptor));
}

static void quad_opposite_set(solver *solver, quad_descriptor desc,
                              vertex_set *set) {
  vertex_set_init_empty(set);

  if (desc == EMPTY_QUAD_DESCRIPTOR)
    return;

  hex_index hex_id = quad_descriptor_hex_index(desc);
  if (hex_id == NO_HEXAHEDRON)
    return;

  quad_index quad_id = quad_descriptor_quad_index(desc);

  vertex_index v = solver->quads[4*quad_id + 0];
  vertex_index u = solver->quads[4*quad_id + 2];
  vertex_index w = solver->quads[4*quad_id + 1];
  vertex_index z = solver->quads[4*quad_id + 3];

  for (hex_vertex_index i = 0; i < 8; i++) {
    vertex_index other = solver->hexes[8*hex_id + i];
    if (other != v && other != u && other != w && other != z)
      vertex_set_insert(set, other);
  }
}

static bool quad_has_same_orientation(vertex_index *q_a, vertex_index *q_b) {
  quad_vertex_index i0;
  for (i0 = 0; q_a[i0] != q_b[0]; i0++);
  return q_b[1] == q_a[(i0 + 1) % 4];
}

static bool finish_quad(solver *solver, fixed_point_delta *delta,
                        hex_index hex_id,
                        vertex_index a, vertex_index b,
                        vertex_index c, vertex_index d) {
  quad_descriptor quad = solver->quad_descriptors[
    diagonal_index(solver->num_vertices, a, b)];
  quad_descriptor other_quad = solver->quad_descriptors[
    diagonal_index(solver->num_vertices, c, d)];

  if (((quad == EMPTY_QUAD_DESCRIPTOR) ^
       (other_quad == EMPTY_QUAD_DESCRIPTOR)) ||
      quad_descriptor_first_vertex(quad) == VERTEX_SET_MAX ||
      quad_descriptor_first_vertex(other_quad) == VERTEX_SET_MAX ||
      (quad != EMPTY_QUAD_DESCRIPTOR &&
       other_quad != EMPTY_QUAD_DESCRIPTOR &&
       quad_descriptor_quad_index(quad) !=
       quad_descriptor_quad_index(other_quad)))
    return false;

  vertex_index vertices[4] = {a, d, b, c};

  if (quad != EMPTY_QUAD_DESCRIPTOR &&
      quad_has_same_orientation(vertices, solver->quads +
                                4 * quad_descriptor_quad_index(quad))) {
    return false;
  }

  quad_index id;
  if (quad == EMPTY_QUAD_DESCRIPTOR) {
    if (solver->max_quad_id != solver->num_quads) {
      delta->max_quad_change++;
      id = solver->max_quad_id++;
      quad = make_quad_descriptor(id, hex_id, c, d);
      solver->quad_descriptors[
        diagonal_index(solver->num_vertices, a, b)] = quad;

      vertex_set_insert(solver->known_quad_diagonals + a, b);
      vertex_set_insert(solver->known_quad_diagonals + b, a);

      quad = make_quad_descriptor(id, hex_id, a, b);
      solver->quad_descriptors[
        diagonal_index(solver->num_vertices, c, d)] = quad;

      vertex_set_insert(solver->known_quad_diagonals + c, d);
      vertex_set_insert(solver->known_quad_diagonals + d, c);

      solver->quads[4*id + 0] = a;
      solver->quads[4*id + 1] = d;
      solver->quads[4*id + 2] = b;
      solver->quads[4*id + 3] = c;

      solver->quad_queue[solver->quad_queue_size++] = id;
    }
    else
      return false;
  }
  else {
    id = quad_descriptor_quad_index(quad);
    vertex_set opposite;
    quad_opposite_set(solver, quad, &opposite);

    for (quad_vertex_index i = 0; i < 4; i++) {
      vertex_set to_keep = solver->domains[i];
      vertex_set_subtract(&to_keep, &opposite);

      vertex_set removed = solver->domains[i];
      vertex_set_subtract(&removed, &to_keep);

      vertex_set_union(&delta->domain_changes[i], &removed);

      solver->domains[i] = to_keep;

      if (vertex_set_empty(&to_keep))
        return false;
    }
  }

  if (solver->num_occurrences[id] == 0)
    return false;

  solver->num_occurrences[id]--;
  delta->occurred[delta->num_occurrences_changes++] = id;

  return true;
}
